Data Manipulation

## read in data and only consider complete data 
## this drops 327 individuals, but BKMR does not handle missing data
nhanes = na.omit(read_csv("./Data/studypop.csv"))
## Parsed with column specification:
## cols(
##   .default = col_double()
## )
## See spec(...) for full column specifications.
## center/scale continous covariates and create indicators for categorical covariates
nhanes$age_z = scale(nhanes$age_cent) ## center and scale age
nhanes$agez_sq = nhanes$age_z^2 ## square this age variable
nhanes$bmicat2 = as.numeric(nhanes$bmi_cat3 == 2) ## 25 <= BMI < 30
nhanes$bmicat3 = as.numeric(nhanes$bmi_cat3 == 3) ## BMI >= 30 (BMI < 25 is the reference)
nhanes$educat1 = as.numeric(nhanes$edu_cat == 1) ## no high school diploma
nhanes$educat3 = as.numeric(nhanes$edu_cat == 3) ## some college or AA degree
nhanes$educat4 = as.numeric(nhanes$edu_cat == 4) ## college grad or above (reference is high schol grad/GED or equivalent)
nhanes$otherhispanic = as.numeric(nhanes$race_cat == 1) ## other Hispanic or other race - including multi-racial
nhanes$mexamerican = as.numeric(nhanes$race_cat == 2) ## Mexican American 
nhanes$black = as.numeric(nhanes$race_cat == 3) ## non-Hispanic Black (non-Hispanic White as reference group)
nhanes$wbcc_z = scale(nhanes$LBXWBCSI)
nhanes$lymphocytes_z = scale(nhanes$LBXLYPCT)
nhanes$monocytes_z = scale(nhanes$LBXMOPCT)
nhanes$neutrophils_z = scale(nhanes$LBXNEPCT)
nhanes$eosinophils_z = scale(nhanes$LBXEOPCT)
nhanes$basophils_z = scale(nhanes$LBXBAPCT)
nhanes$lncotinine_z = scale(nhanes$ln_lbxcot) ## to access smoking status, scaled ln cotinine levels


## our y variable - ln transformed and scaled mean telomere length
lnLTL_z = scale(log(nhanes$TELOMEAN)) 

## our Z matrix
mixture = with(nhanes, cbind(LBX074LA, LBX099LA, LBX118LA, LBX138LA, LBX153LA, LBX170LA, LBX180LA,
                             LBX187LA, LBX194LA, LBXHXCLA, LBXPCBLA, LBXD03LA, LBXD05LA, LBXD07LA,
                             LBXF03LA, LBXF04LA, LBXF05LA, LBXF08LA)) 
lnmixture = apply(mixture, 2, log)
lnmixture_z = scale(lnmixture)
colnames(lnmixture_z) = c(paste0("PCB",c(74, 99, 118, 138, 153, 170, 180, 187, 194, 169, 126)), 
                           paste0("Dioxin",1:3), paste0("Furan",1:4)) 

## our X matrix
covariates = with(nhanes, cbind(age_z, agez_sq, male, bmicat2, bmicat3, educat1, educat3, educat4, 
                                 otherhispanic, mexamerican, black, wbcc_z, lymphocytes_z, monocytes_z, 
                                 neutrophils_z, eosinophils_z, basophils_z, lncotinine_z)) 


### create knots matrix for Gaussian predictive process (to speed up BKMR with large datasets)
set.seed(10)
knots100     <- fields::cover.design(lnmixture_z, nd = 100)$design
#save(knots100, knots100.PCB, file = "./Supervised/BKMR/saved_model/NHANES_knots100.RData") 

Fit Models

Group VS Fit with all Exposures using GPP and 100 Knots

### Group VS fit with all exposures using GPP and 100 knots 
load("./Supervised/BKMR/saved_model/NHANES_knots100.RData")

##### fit BKMR models WITH Gaussian predictive process using 100 knots

### Group VS fit with all exposures using GPP and 100 knots 
set.seed(1000)
## Note: the commented out statement is what generated the saved model fits. 
## For this lab we are only going to generate 100 MCMC samples to get a sense
##   for what the program outputs.  

#fit_gvs_knots100 <-  kmbayes(y=lnLTL_z, Z=lnmixture_z, X=covariates, iter=100000, verbose=TRUE, varsel=TRUE, 
#                             groups=c(rep(1,times=2), 2, rep(1,times=6), rep(3,times=2),rep(2,times=7)), knots=knots100)

temp <-  kmbayes(y=lnLTL_z, Z=lnmixture_z, X=covariates, iter=100, verbose=TRUE, varsel=TRUE, 
                             groups=c(rep(1,times=2), 2, rep(1,times=6), rep(3,times=2),
                                      rep(2,times=7)), knots=knots100)
## Iteration: 10 (10% completed; 2.29639 secs elapsed)
## Acceptance rates for Metropolis-Hastings algorithm:
##               param      rate
## 1            lambda 0.2222222
## 2 r/delta (overall) 0.7777778
## 3 r/delta  (move 1) 1.0000000
## 4 r/delta  (move 2) 0.6000000
## 5 r/delta  (move 3)       NaN
## Iteration: 20 (20% completed; 4.26684 secs elapsed)
## Acceptance rates for Metropolis-Hastings algorithm:
##               param      rate
## 1            lambda 0.1578947
## 2 r/delta (overall) 0.4736842
## 3 r/delta  (move 1) 0.5555556
## 4 r/delta  (move 2) 0.3750000
## 5 r/delta  (move 3) 0.5000000
## Iteration: 30 (30% completed; 6.21747 secs elapsed)
## Acceptance rates for Metropolis-Hastings algorithm:
##               param      rate
## 1            lambda 0.1034483
## 2 r/delta (overall) 0.4137931
## 3 r/delta  (move 1) 0.4166667
## 4 r/delta  (move 2) 0.4166667
## 5 r/delta  (move 3) 0.4000000
## Iteration: 40 (40% completed; 8.30155 secs elapsed)
## Acceptance rates for Metropolis-Hastings algorithm:
##               param       rate
## 1            lambda 0.07692308
## 2 r/delta (overall) 0.46153846
## 3 r/delta  (move 1) 0.42857143
## 4 r/delta  (move 2) 0.47368421
## 5 r/delta  (move 3) 0.50000000
## Iteration: 50 (50% completed; 10.31038 secs elapsed)
## Acceptance rates for Metropolis-Hastings algorithm:
##               param       rate
## 1            lambda 0.06122449
## 2 r/delta (overall) 0.38775510
## 3 r/delta  (move 1) 0.46666667
## 4 r/delta  (move 2) 0.36000000
## 5 r/delta  (move 3) 0.33333333
## Iteration: 60 (60% completed; 12.33724 secs elapsed)
## Acceptance rates for Metropolis-Hastings algorithm:
##               param       rate
## 1            lambda 0.05084746
## 2 r/delta (overall) 0.40677966
## 3 r/delta  (move 1) 0.43750000
## 4 r/delta  (move 2) 0.39393939
## 5 r/delta  (move 3) 0.40000000
## Iteration: 70 (70% completed; 14.61572 secs elapsed)
## Acceptance rates for Metropolis-Hastings algorithm:
##               param      rate
## 1            lambda 0.1304348
## 2 r/delta (overall) 0.4057971
## 3 r/delta  (move 1) 0.4500000
## 4 r/delta  (move 2) 0.4166667
## 5 r/delta  (move 3) 0.3076923
## Iteration: 80 (80% completed; 16.61848 secs elapsed)
## Acceptance rates for Metropolis-Hastings algorithm:
##               param      rate
## 1            lambda 0.1265823
## 2 r/delta (overall) 0.3924051
## 3 r/delta  (move 1) 0.4545455
## 4 r/delta  (move 2) 0.4047619
## 5 r/delta  (move 3) 0.2666667
## Iteration: 90 (90% completed; 18.54577 secs elapsed)
## Acceptance rates for Metropolis-Hastings algorithm:
##               param      rate
## 1            lambda 0.1460674
## 2 r/delta (overall) 0.3820225
## 3 r/delta  (move 1) 0.4000000
## 4 r/delta  (move 2) 0.4000000
## 5 r/delta  (move 3) 0.3157895
## Iteration: 100 (100% completed; 20.49031 secs elapsed)
## Acceptance rates for Metropolis-Hastings algorithm:
##               param      rate
## 1            lambda 0.1414141
## 2 r/delta (overall) 0.3535354
## 3 r/delta  (move 1) 0.3571429
## 4 r/delta  (move 2) 0.3673469
## 5 r/delta  (move 3) 0.3181818
## The following statement saved the model fits using 100,000 MCMC samples. We won't save
##    our fit based on 100 samples. Rather we will load in the results from 100,000. 

#save(fit_gvs_knots100,file="bkmr_NHANES_gvs_knots100.RData")
load("./Supervised/BKMR/saved_model/bkmr_NHANES_gvs_knots100.RData")
summary(fit_gvs_knots100)
## Fitted object of class 'bkmrfit'
## Iterations: 1e+05 
## Outcome family: gaussian  
## Model fit on: 2018-12-11 19:32:59 
## Running time:  6.46064 hours 
## 
## Acceptance rates for Metropolis-Hastings algorithm:
##               param       rate
## 1            lambda 0.08499085
## 2 r/delta (overall) 0.24789248
## 3 r/delta  (move 1) 0.33862879
## 4 r/delta  (move 2) 0.24140731
## 5 r/delta  (move 3) 0.16400060
## 
## Parameter estimates (based on iterations 50001-100000):
##        param     mean      sd    q_2.5   q_97.5
## 1      beta1 -0.53023 0.04356 -0.61920 -0.44620
## 2      beta2 -0.01427 0.03218 -0.07732  0.04899
## 3      beta3 -0.15517 0.06135 -0.27624 -0.03525
## 4      beta4 -0.01796 0.07069 -0.15548  0.12091
## 5      beta5 -0.06813 0.07909 -0.22397  0.08561
## 6      beta6 -0.08602 0.07921 -0.24138  0.06977
## 7      beta7  0.06920 0.07919 -0.08595  0.22516
## 8      beta8 -0.01479 0.08959 -0.19039  0.16028
## 9      beta9  0.17799 0.10938 -0.03762  0.39132
## 10    beta10  0.04093 0.08263 -0.12039  0.20275
## 11    beta11  0.19093 0.08279  0.02932  0.35303
## 12    beta12 -0.05792 0.03329 -0.12330  0.00741
## 13    beta13 -1.10570 3.06012 -7.08628  4.92819
## 14    beta14 -0.31742 0.74194 -1.76554  1.14300
## 15    beta15 -1.18773 3.39199 -7.81900  5.51197
## 16    beta16 -0.27092 0.78130 -1.80130  1.26774
## 17    beta17 -0.02009 0.16103 -0.33394  0.29761
## 18    beta18  0.05246 0.03224 -0.01091  0.11560
## 19 sigsq.eps  0.74942 0.03396  0.68556  0.81925
## 20        r1  0.00087 0.00688  0.00000  0.01275
## 21        r2  0.00136 0.02158  0.00000  0.01467
## 22        r3  0.00124 0.00812  0.00000  0.01664
## 23        r4  0.00109 0.00624  0.00000  0.01575
## 24        r5  0.00131 0.00832  0.00000  0.01555
## 25        r6  0.00202 0.01103  0.00000  0.02380
## 26        r7  0.00153 0.00783  0.00000  0.01966
## 27        r8  0.00101 0.00778  0.00000  0.01333
## 28        r9  0.00095 0.00613  0.00000  0.01454
## 29       r10  0.00652 0.02102  0.00000  0.05517
## 30       r11  0.00783 0.01333  0.00000  0.04042
## 31       r12  0.00009 0.00166  0.00000  0.00000
## 32       r13  0.00021 0.00481  0.00000  0.00000
## 33       r14  0.00015 0.00230  0.00000  0.00000
## 34       r15  0.01832 0.02523  0.00000  0.07990
## 35       r16  0.00018 0.00230  0.00000  0.00000
## 36       r17  0.00028 0.00350  0.00000  0.00000
## 37       r18  0.00046 0.00590  0.00000  0.00000
## 38    lambda  2.14163 3.51586  0.14823 11.74103
## 
## Posterior inclusion probabilities:
##    variable group groupPIP condPIP
## 1     PCB74     1  0.41952 0.09320
## 2     PCB99     1  0.41952 0.12710
## 3    PCB118     2  0.86570 0.05656
## 4    PCB138     1  0.41952 0.11709
## 5    PCB153     1  0.41952 0.13816
## 6    PCB170     1  0.41952 0.17263
## 7    PCB180     1  0.41952 0.15098
## 8    PCB187     1  0.41952 0.09926
## 9    PCB194     1  0.41952 0.10159
## 10   PCB169     3  0.61998 0.35298
## 11   PCB126     3  0.61998 0.64702
## 12  Dioxin1     2  0.86570 0.00534
## 13  Dioxin2     2  0.86570 0.00880
## 14  Dioxin3     2  0.86570 0.00737
## 15   Furan1     2  0.86570 0.88037
## 16   Furan2     2  0.86570 0.01047
## 17   Furan3     2  0.86570 0.01243
## 18   Furan4     2  0.86570 0.01867

Posterior Inclusion Probabilities (PIPs)

## obtain posterior inclusion probabilities (PIPs)
ExtractPIPs(fit_gvs_knots100)
##    variable group groupPIP     condPIP
## 1     PCB74     1  0.41952 0.093201754
## 2     PCB99     1  0.41952 0.127097635
## 3    PCB118     2  0.86570 0.056555389
## 4    PCB138     1  0.41952 0.117086194
## 5    PCB153     1  0.41952 0.138157895
## 6    PCB170     1  0.41952 0.172625858
## 7    PCB180     1  0.41952 0.150982075
## 8    PCB187     1  0.41952 0.099256293
## 9    PCB194     1  0.41952 0.101592296
## 10   PCB169     3  0.61998 0.352979128
## 11   PCB126     3  0.61998 0.647020872
## 12  Dioxin1     2  0.86570 0.005336722
## 13  Dioxin2     2  0.86570 0.008802125
## 14  Dioxin3     2  0.86570 0.007369759
## 15   Furan1     2  0.86570 0.880374264
## 16   Furan2     2  0.86570 0.010465519
## 17   Furan3     2  0.86570 0.012429248
## 18   Furan4     2  0.86570 0.018666975

Data Visualization

Correlation Matrix

### correlation matrix
cor.Z = cor(lnmixture_z, use = "complete.obs")

#pdf(file= "./Supervised/BKMR/saved_model/cor_nhanes.pdf", width=12, height=12)
corrplot.mixed(cor.Z, upper = "ellipse", lower.col = "black")

#dev.off()

Model

Trace Plots

### change this for each model you fit and then rerun the code from here to the bottom
modeltoplot      <- fit_gvs_knots100   ## name of model object
modeltoplot.name <- "fit_gvs_knots100" ## name of model for saving purposes
plot.name        <- "gvs_knots100"     ## part that changed in plot name 
Z                <- lnmixture_z        ## Z matrix to match what was used in model

### values to keep after burnin/thin
sel<-seq(50001,100000,by=50)


### access convergence with traceplots 
TracePlot(fit = modeltoplot, par = "beta", sel=sel)

TracePlot(fit = modeltoplot, par = "sigsq.eps", sel=sel)

par(mfrow=c(2,2))
TracePlot(fit = modeltoplot, par = "r", comp = 1, sel=sel)
TracePlot(fit = modeltoplot, par = "r", comp = 2, sel=sel)
TracePlot(fit = modeltoplot, par = "r", comp = 3, sel=sel)
TracePlot(fit = modeltoplot, par = "r", comp = 4, sel=sel)

par(mfrow=c(1,1))

Univariable Exposure-Response Functions

#### create dataframes for ggplot (this takes a little while to run)

## We will not run these 6 commands that generate all the data for the summaries of
## the exposure-response function. Instead we will read in the previously generated results 
## just to save a little time. 

#pred.resp.univar <- PredictorResponseUnivar(fit = modeltoplot, sel=sel, method="approx")
#pred.resp.bivar  <- PredictorResponseBivar(fit = modeltoplot,  min.plot.dist = 1, sel=sel, method="approx")
#pred.resp.bivar.levels <- PredictorResponseBivarLevels(pred.resp.df = pred.resp.bivar, Z = Z,
#                                                          both_pairs = TRUE, qs = c(0.25, 0.5, 0.75))
#risks.overall <- OverallRiskSummaries(fit = modeltoplot, qs = seq(0.25, 0.75, by = 0.05), q.fixed = 0.5, method = "approx",sel=sel)
#risks.singvar <- SingVarRiskSummaries(fit = modeltoplot, qs.diff = c(0.25, 0.75),
#                                        q.fixed = c(0.25, 0.50, 0.75), method = "approx")
#risks.int <- SingVarIntSummaries(fit = modeltoplot, qs.diff = c(0.25, 0.75), qs.fixed = c(0.25, 0.75))

#Save the data need to plot the summaries of interest generated 
#save(pred.resp.univar, pred.resp.bivar, pred.resp.bivar.levels, risks.overall, risks.singvar, risks.int, file=paste0("./Supervised/BKMR/saved_model/", modeltoplot.name,"_plots.RData"))

# Load in the results, which were computed previously and saved using the command just above this

load(paste0("./Supervised/BKMR/saved_model/", modeltoplot.name,"_plots.RData"))
pred.resp.univar %>% 
  mutate(variable = fct_recode(variable, "PCB 74" = "PCB74",
                                "PCB 99" = "PCB99",
                                "PCB 118" = "PCB118",
                                "PCB 138" = "PCB138",
                                "PCB 153" = "PCB153",
                                "PCB 170" = "PCB170",
                                "PCB 180" = "PCB180",
                                "PCB 187" = "PCB187",
                                "PCB 194" = "PCB194",
                                "1,2,3,6,7,8-hxcdd" = "Dioxin1",
                                "1,2,3,4,6,7,8-hpcdd" = "Dioxin2",
                               "1,2,3,4,6,7,8,9-ocdd" =  "Dioxin3",
                               "2,3,4,7,8-pncdf" =  "Furan1",
                               "1,2,3,4,7,8-hxcdf" =  "Furan2",
                               "1,2,3,6,7,8-hxcdf" =  "Furan3",
                               "1,2,3,4,6,7,8-hxcdf" =  "Furan4",
                               "PCB 169" =  "PCB169",
                                "PCB 126" = "PCB126")) %>% 
  ggplot(aes(z, est, ymin = est - 1.96*se, ymax = est + 1.96*se)) + 
  geom_smooth(stat = "identity") + labs(y = "Estimate", x = "Exposure") + 
  facet_wrap(~ variable) + theme_bw() +
  theme(strip.background = element_rect(fill = "white"))

Bivariable Exposure-Response Functions

pred.resp.bivar %>% 
  mutate(variable1 = fct_recode(variable1, "PCB 74" = "PCB74",
                                "PCB 99" = "PCB99",
                                "PCB 118" = "PCB118",
                                "PCB 138" = "PCB138",
                                "PCB 153" = "PCB153",
                                "PCB 170" = "PCB170",
                                "PCB 180" = "PCB180",
                                "PCB 187" = "PCB187",
                                "PCB 194" = "PCB194",
                                "1,2,3,6,7,8- hxcdd" = "Dioxin1",
                                "1,2,3,4,6,7,8- hpcdd" = "Dioxin2",
                               "1,2,3,4,6,7,8,9- ocdd" =  "Dioxin3",
                               "2,3,4,7,8- pncdf" =  "Furan1",
                               "1,2,3,4,7,8- hxcdf" =  "Furan2",
                               "1,2,3,6,7,8- hxcdf" =  "Furan3",
                               "1,2,3,4,6,7,8- hxcdf" =  "Furan4",
                               "PCB 169" =  "PCB169",
                                "PCB 126" = "PCB126"),
       variable2 = fct_recode(variable2, "PCB 74" = "PCB74",
                                "PCB 99" = "PCB99",
                                "PCB 118" = "PCB118",
                                "PCB 138" = "PCB138",
                                "PCB 153" = "PCB153",
                                "PCB 170" = "PCB170",
                                "PCB 180" = "PCB180",
                                "PCB 187" = "PCB187",
                                "PCB 194" = "PCB194",
                                "1,2,3,6,7,8- hxcdd" = "Dioxin1",
                                "1,2,3,4,6,7,8- hpcdd" = "Dioxin2",
                               "1,2,3,4,6,7,8,9- ocdd" =  "Dioxin3",
                               "2,3,4,7,8- pncdf" =  "Furan1",
                               "1,2,3,4,7,8- hxcdf" =  "Furan2",
                               "1,2,3,6,7,8- hxcdf" =  "Furan3",
                               "1,2,3,4,6,7,8- hxcdf" =  "Furan4",
                               "PCB 169" =  "PCB169",
                                "PCB 126" = "PCB126")) %>% 
ggplot(aes(z1, z2, fill = est)) + 
  geom_raster() + 
  facet_grid(variable2 ~ variable1, scales = "free", space = "free",
             labeller = labeller(variable1 = label_wrap_gen(5),
                                 variable2 = label_wrap_gen(5))) +
  scale_fill_gradientn(colours=c("#0000FFFF","#FFFFFFFF","#FF0000FF")) +
  xlab("Exposure 1") +
  ylab("Exposure 2") + labs(fill = "Estimate") +
  ggtitle("h(Exposure 1, Exposure 2)") + theme_bw() +
  theme(strip.background = element_rect(fill = "white"),
        panel.spacing = unit(0.05, "lines"),
        legend.position = "bottom")

pred.resp.bivar.levels %>% 
mutate(variable1 = fct_recode(variable1, "PCB 74" = "PCB74",
                                "PCB 99" = "PCB99",
                                "PCB 118" = "PCB118",
                                "PCB 138" = "PCB138",
                                "PCB 153" = "PCB153",
                                "PCB 170" = "PCB170",
                                "PCB 180" = "PCB180",
                                "PCB 187" = "PCB187",
                                "PCB 194" = "PCB194",
                                "1,2,3,6,7,8- hxcdd" = "Dioxin1",
                                "1,2,3,4,6,7,8- hpcdd" = "Dioxin2",
                               "1,2,3,4,6,7,8,9- ocdd" =  "Dioxin3",
                               "2,3,4,7,8- pncdf" =  "Furan1",
                               "1,2,3,4,7,8- hxcdf" =  "Furan2",
                               "1,2,3,6,7,8- hxcdf" =  "Furan3",
                               "1,2,3,4,6,7,8- hxcdf" =  "Furan4",
                               "PCB 169" =  "PCB169",
                                "PCB 126" = "PCB126"),
       variable2 = fct_recode(variable2, "PCB 74" = "PCB74",
                                "PCB 99" = "PCB99",
                                "PCB 118" = "PCB118",
                                "PCB 138" = "PCB138",
                                "PCB 153" = "PCB153",
                                "PCB 170" = "PCB170",
                                "PCB 180" = "PCB180",
                                "PCB 187" = "PCB187",
                                "PCB 194" = "PCB194",
                                "1,2,3,6,7,8- hxcdd" = "Dioxin1",
                                "1,2,3,4,6,7,8- hpcdd" = "Dioxin2",
                               "1,2,3,4,6,7,8,9- ocdd" =  "Dioxin3",
                               "2,3,4,7,8- pncdf" =  "Furan1",
                               "1,2,3,4,7,8- hxcdf" =  "Furan2",
                               "1,2,3,6,7,8- hxcdf" =  "Furan3",
                               "1,2,3,4,6,7,8- hxcdf" =  "Furan4",
                               "PCB 169" =  "PCB169",
                                "PCB 126" = "PCB126")) %>% 
  ggplot(aes(z1, est)) + 
  geom_smooth(aes(col = quantile), stat = "identity") + 
  facet_grid(variable2 ~ variable1, scales = "free", space = "free",
             labeller = labeller(variable1 = label_wrap_gen(5),
                                 variable2 = label_wrap_gen(5))) +
  ggtitle("h(Exposure 1 | Quantiles of Exposure 2)") +
  xlab("Exposure 1") + theme_bw() + labs(col = "Quantile", y = "Estimate") +
  theme(strip.background = element_rect(fill = "white"),
        panel.spacing = unit(0.05, "lines"),
        legend.position = "bottom")
## Warning: Removed 2190 rows containing missing values (geom_smooth).

Overall Mixture Effect

ggplot(risks.overall, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) +  
  geom_hline(yintercept=00, linetype="dashed", color="gray") + 
  geom_pointrange() + theme_bw() +
  labs(x = "Quantile", y = "Estimate")

Single Variable

risks.singvar %>% mutate(variable = fct_recode(variable, "PCB 74" = "PCB74",
                                "PCB 99" = "PCB99",
                                "PCB 118" = "PCB118",
                                "PCB 138" = "PCB138",
                                "PCB 153" = "PCB153",
                                "PCB 170" = "PCB170",
                                "PCB 180" = "PCB180",
                                "PCB 187" = "PCB187",
                                "PCB 194" = "PCB194",
                                "1,2,3,6,7,8-hxcdd" = "Dioxin1",
                                "1,2,3,4,6,7,8-hpcdd" = "Dioxin2",
                               "1,2,3,4,6,7,8,9-ocdd" =  "Dioxin3",
                               "2,3,4,7,8-pncdf" =  "Furan1",
                               "1,2,3,4,7,8-hxcdf" =  "Furan2",
                               "1,2,3,6,7,8-hxcdf" =  "Furan3",
                               "1,2,3,4,6,7,8-hxcdf" =  "Furan4",
                               "PCB 169" =  "PCB169",
                                "PCB 126" = "PCB126")) %>% 
ggplot(aes(variable, est, ymin = est - 1.96*sd,  ymax = est + 1.96*sd, col = q.fixed)) + 
  geom_hline(aes(yintercept=0), linetype="dashed", color="gray") +  
  geom_pointrange(position = position_dodge(width = 0.75)) +  coord_flip() + 
  theme_bw() +
  labs(x = "", y = "Estimate", col = "Fixed Quantile")

Single Variable Interaction Terms

risks.int %>% mutate(variable = fct_recode(variable, "PCB 74" = "PCB74",
                                "PCB 99" = "PCB99",
                                "PCB 118" = "PCB118",
                                "PCB 138" = "PCB138",
                                "PCB 153" = "PCB153",
                                "PCB 170" = "PCB170",
                                "PCB 180" = "PCB180",
                                "PCB 187" = "PCB187",
                                "PCB 194" = "PCB194",
                                "1,2,3,6,7,8-hxcdd" = "Dioxin1",
                                "1,2,3,4,6,7,8-hpcdd" = "Dioxin2",
                               "1,2,3,4,6,7,8,9-ocdd" =  "Dioxin3",
                               "2,3,4,7,8-pncdf" =  "Furan1",
                               "1,2,3,4,7,8-hxcdf" =  "Furan2",
                               "1,2,3,6,7,8-hxcdf" =  "Furan3",
                               "1,2,3,4,6,7,8-hxcdf" =  "Furan4",
                               "PCB 169" =  "PCB169",
                                "PCB 126" = "PCB126")) %>% 
ggplot(aes(variable, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) + 
  geom_pointrange(position = position_dodge(width = 0.75)) + 
  geom_hline(yintercept = 0, lty = 2, col = "gray") + coord_flip() + theme_bw() +
  labs(x = "", y = "Estimate")